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1. Introduction 

QCD at finite temperature (T) and/or chemical potential (/x) is of fundamental importance, 
since it describes relevant features of particle physics in the early universe, in neutron stars 
and in heavy ion collisions (for a clear introduction see [1]). QCD is asymptotically free, 
thus its high T and high density phases are dominated by partons (quarks and gluons) 
as degrees of freedom rather than hadrons. In this quark-gluon plasma (QGP) phase the 
symmetries of QCD are restored. In addition, recently a particularly interesting, rich phase 
structure has been conjectured for QCD at finite T and /x [2-5]. 

Much effort has been devoted recently to heavy ion collisions at CERN and Brookhaven 
in order to experimentally detect the quark-gluon plasma. Clearly, a theoretical under- 
standing of the underlying physics is of extreme importance. Extensive lattice QCD calcu- 
lations were carried out to give first principle answers e.g. for the transition^ temperature 
(Tc) and for the equation of state (EoS). For recent reviews see [6-14]. Unfortunately, all 
these results are at vanishing chemical potential, whereas the experiments are carried out 
at non- vanishing /i values. 

Thus, the main goal of the present paper is to determine the EoS at finite temperature 
and chemical potential. 

^We use the expression "transition" if we do not want to specify whether we deal with a phase transition 
or a crossover. 
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QCD at finite ^ can be formulated on the lattice [15, 16]; however, standard Monte- 
Carlo techniques cannot be used at 7^ 0. The reason is that for non-vanishing real /j, the 
functional measure - thus, the determinant of the Euclidean Dirac operator - is complex. 
This fact spoils any Monte-Carlo technique based on importance sampling. Several pro- 
posals were studied to solve the problem. Unfortunately, none of them was able to give 
the EoS at non- vanishing ^. 

In a recent paper two of us proposed a new method, the so-called overlap improv- 
ing multi-parameter rcwcighting technique [17], to study lattice QCD and give the phase 
boundary at finite T and /x. The idea was to produce an ensemble of QCD configurations 
at /Lt=0 and at T = Tc. Then the Boltzmann weights of these configurations at 7^ 
and at T lowered to the transition temperatures were determined at this non-vanishing fj, 
using [18,19]. Since transition configurations were reweighted to transition configurations 
a much better overlap was observed than by reweighting pure hadronic configurations to 
transition ones [20]. We also emphasized that for small n the technique works for tem- 
peratures both below and above the transition temperature. We generalized the overlap 
improving multi-parameter reweighting method to arbitrary number of staggered quarks 
and applied it to the Uf = 2 + 1 case [21]. Based on the volume (V) dependence of the 
Lee- Yang zeros of the partition function we determined the endpoint ^ of QCD with semi- 
realistic masses on A^^ = 4 lattices. We obtained ~ 160 McV for the endpoint temperature 
and !=s! 700 MeV for the endpoint baryonic chemical potential. Note that using a Taylor 
expansion around /U = 0, T 7^ for small chemical potentials could be also seen as a vari- 
ant of the multi-parameter reweighting method, which can be used to determine hadron 
masses [24], thermal properties [25] and even to obtain the EoS as was done in [26] for 
two flavours. The Taylor expansion method was also used to evaluate the pressure on 
quenched configurations [27]. Recently, simulations at imaginary chemical potentials and 
analytic continuation were also used to determine the phase boundary on the jJL — T plane 
for 2,3 and 4 flavour staggered QCD [28-30]. 

In this paper we suggest a technique by which the EoS is determined on the line of 
constant physics (LCP) ^. An LCP can be defined ^ by a fixed ratio of the strange quark 
mass {nis) and light quark masses (mud) to the /U=0 transition temperature (Tc). The /i = 
LCP results are compared with earlier "non-LCP techniques" ( [31-33]) , which calculate 
the EoS at fixed m^a (in this approach the increase of the temperature - thus the decrease 
of the lattice spacings ''a" - results in an increase of the quark mass). We comment on 
the differences. Our parameter choice corresponds approximately to the physical strange 
quark mass. However, the ratio of the pion mass (m^r) and the mass of the rho meson 
(mp) is around 0.5 — 0.75, which is roughly 3 times larger than its physical value. The 
temperature dependence of the EoS is studied in a wide range. In our lattice analysis 

^The same combination of the multi-parameter reweighting and the Lee- Yang technique was successfully 
used on quite large lattices (upto spatial volumes of 60^) to locate the endpoint of the electroweak phase 
transition [22,23]. 

^The importance of using LCP's is explained in Section 3 below. 

''Our first choice for an LCP is given by the bare masses, later we determine the LCP using renormalized 
quantites. 
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we use 2 + 1 flavour QCD with dynamical staggered quarks. The determination of the 
equation of state at finite chemical potential requires several observables (e.g. expectation 
values of the plaquette {PI), or the chiral condensate {^^)) at non- vanishing // values. In 
order to obtain these quantities we use the overlap improving multi-parameter reweighting 
technique of Ref. [17]. We employ the integral method [34] to calculate the pressure. The 
energy density (e) can be obtained by combining the results on the pressure and on the 
"interaction measure" (e — 3p). Using the pressure, one can directly determine the quark 
number density (n). 

The paper is organized as follows. In Section 2 we summarize the lattice parameters. In 
Section 3 we present the technique by which the lines of constant physics can be determined. 
Section 4 presents the equation of state at vanishing chemical potential. Sections 5 and 
6 deal with the question how to reweight into the region of 7^ and how to estimate 
the error of the reweightcd quantities. In Section 7 we give the equation of state for non- 
vanishing chemical potential and temperature. Those who are not interested in the details 
of the lattice techniques should simply omit Sections 2-6 and jump to Section 7, or refer 
to [35] or [36]. Finally, Section 8 contains a summary and the conclusions. 

2. Lattice parameters 

In this paper we use 2 + 1 flavour dynamical QCD with unimproved staggered action. We 

study the system at several gauge couplings and masses. Simulations arc done for the equa- 
tion of state along two different lines of constant physics and at 14 different temperatures. 
Following the experiences of previous studies on the EoS the simulation points are dis- 
tributed in a way that they are denser around the transition temperature than elsewhere. 
The temperature range spans up to 3Tc. In physical units our parameters correspond to 
pion to rho mass ratio of rriT^Imp 0.5 — 0.75 and lattice spacings of a 0.12 — 0.35 fm. 
See the next section for more details on the measured quantities (masses, string tension: 
^/a and potential scale: Rq,Ri) along the lines of constant physics. 

The finite temperature contributions to the EoS are obtained on 4- 8^, 4- 10-^ and 4-12''^ 
lattices, which can be used to extrapolate into the thermodynamical limit (we usually call 
them hot lattices) . On these lattices we determine not only the usual observables (plaque- 
tte, Polyakov line, chiral condensates) but also the determinant of the fermion matrix and 
the baryon density (n^) at finite n. 10000 — 20000 trajectories are simulated at each bare 
parameter set. Plaqucttcs, Polyakov lines and the chiral condensates are measured at each 
trajectory whereas the CPU demanding determinants and related quantities arc evaulatcd 
at every 30 trajectories. For our parameters the CPU time used for the production of 
configurations is of the same order of magnitude than the CPU time used for calculating 
the determinants. 

Zero temperature runs are done on 12^ and 14^ • 24 lattices (we usually call them 
cold lattices). 1500 — 2500 trajectories are simulated at each bare parameter set. Hadron 
propagators and Wilson loops are measured at every 10 trajectories. Masses and the static 
potential are determined by correlated fits. The proper fitting interval is chosen to obtain 
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Figure 1: (a) The variation of the plaquette as a function of the step-size squared at /? = 5.4 
and nriqa — 0.064. The boxes come from the cold and the triangles from the hot lattices. The cold 
system has a stronger step-size dependence. As it can be seen step-sizes smaller than niudj'^ result 
in smaller systematic uncertainties than the statistical error, (b) The same for the chiral condensate 
of the light quarks. Their step-size dependence is negligible. 

minimal x^/d-o.f., and/or x^/d.o.f.~l. In order to avoid instability when inverting the 
correlation matrix we use the smeared eigenvalue technique [37]. 

Our simulations employ the standard R-algorithm [38]. The length of one trajectory 
is unity. The step-size dependence of different relevant quantities (expectation values of 
the plaquette (-P/), or chiral condensate (^^)) was already analysed in the literature (see 
e.g. [31,39]). We also studied the systematic uncertainties associated with the step-size 
dependence of the results. Similarly to other groups we found that these effects are much 
more pronounced for cold lattices. The most sensitive quantity is the plaquette difference 
between hot and cold lattices. This difference can be of the order of the step-size error. 
Figure |l| shows the plaquette and the chiral condensate of the light quarks as a function of 
the step-size squared, which is the leading error in the R-algorithm [38]. Based on these 
experiences the step-size is chosen to be m„rf/2. The errors that are introduced this way 
are 0.2 % for zero temperature lattices and 0.1 % for finite temperature ones. 

Statistical uncertainities are determined by jackknife analysis with 20-50 jackknife 
samples. 

Since we usually move along the line of constant physics by changing the lattice spac- 
ing a and keeping the masses fixed we will explicitely write out the lattice spacing a in 
our formulas. In this paper we study lattices with isotropic couplings. In this case the 
lattice spacing is the same in the temporal and spatial directions. We write [Ib for the 
baryonic chemical potential, whereas for the quark chemical potential (n, d quarks) we use 
the notation /i. Similarly, the baryon density is denoted by and the light quark density 
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by n. 

3. Lines of constant physics (LCP) at = 

In this section we discuss the role of LCP when determining the EoS in pure gauge theory 
and in dynamical QCD. After that we determine the lines of constant physics, along which 
our simulations are done. 

The EoS has been determined for pure gauge (quenched) theory on the lattice already 
in the early years of lattice QCD. Continuum extrapolations have been made by using 
several lattice actions. The results are in agreement within errorbars [40-42]. Results 
obtained by the integral method and the derivative method were also compared [43,44]. 

In order to detetermine the temperature T = l/[Nta{(3)] of the pure gauge theory, 
we have to compute the lattice spacing (o) as a function of the gauge coupling (/?). Note 
that when changing the coupling of the pure gauge system one automatically remains 
(upto scaling violations) on the line of constant physics. The situation is quite different in 
full QCD. In the d dimensional space of the bare parameters one defines d appropriately 
chosen quantities. The LCP is given by d — 1 constraints and it is parametrized by a non- 
constrained combination of the above quantities. For the 2 + 1 flavour staggered action we 
have three bare parameters (/?, and two masses, m„rf,ms). Thus, we need two constraints. 
There are several possibilities for these constraints and consequently there are many ways 
to define an LCP. A convenient choice for two of the three quantities can be the bare quark 
masses {mud and m^). A more physical possibility is to use the pion and kaon masses 
(m,r,mi^). There are several options for the third quantity. It can be the rho mass (rrip), 
the string tension {\/^), the Ro,Ri scales of the potential ( [45], and also [46]) or the 
transition temperature. Fixing the ratios of the third quantitiy to the first two gives two 
constraints and the third quantity in lattice units fixes the scale along the LCP. 

In our analysis first we use the bare quark masses {rriud and m<j) and the transition 
temperature to define an LCP. In this paper we use two^ LCP's (LCPi and LCP2). The 
conditions 

m„rf = 0.48re = 0.48/(iVta) and m, = 2.08 • m„rf 
mud = 0.384rc = 0.384/(7Vta) and = 2.08 • m^d ' 

are used as the constraints for LCPi and LCP2, respectively^. For both LCP's we deter- 
mined four different transition couplings (/3c) by susceptibility peaks on Nt = l/{Tca) = 4, 
6, 8 and 10 lattices with spatial extensions Ns > 'iNt/2 and quark masses given by eq. (p.l|). 
The quark masses or the transition gauge couplings can be used to parametrize the LCP's. 

Note that mg/m^d is the same for both of our LCP's. This relationship does not change 
when we interpolate between the two LCP's or perform reweighting in some parameters. 
Thus, for transparency we usually do not specify the quark fiavours when discussing the 

^As we will see later, two LCP's are needed for the determination of the EoS at finite chemical potential. 
®Our light quark masses are heavier than the physical ones, whereas the strange quark mass is at its 
approximate physical value. 
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physical scale. We simply speak about "quark mass parameter" or "quark mass". The 
flavour dependence is indicated explicitely only when necessary. 

It is instructive to define LCP's by using renormalized quantities (LCP*'s, namely 
LCPJ and LCP2) and to test scaling violation. One can calculate mp,m,r) Ro and a on 
T = lattices {V = 14^ • 24) using the bare parameters obtained for LCPi and LCP2. 
By measuring these quantities at somewhat different quark masses gives the possibility to 
define by linear interpolations the LCP*'s. Our results for LCP2 and LCP2 are summarized 
in Table |l] and Table |2[ The definition of LCP2 was m-n^/nip ~ 0.626. As it can be seen 
different definitions of the scales deviate from each other only by ~ 8 %. Table [l| and Table 
|2| contains the T = results for the "non-LCP approach", too (see next section). 



Nt 
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ma 
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Ri/a 


nipa 




4 


5.271 


0.096 


0.5889(9) 


2.02(1) 


1.483(2) 


1.421(1) 


0.78514(6) 
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5.4 


0.064 


0.3686(6) 


3.19(5) 


2.34(2) 


1.048(1) 


0.6805(1) 
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5.5 


0.048 


0.2697(5) 


3.96(9) 


3.64(10) 


0.778(4) 


0.5595(13) 


10 


5.58 


0.0384 


0.2358(4) 


4.44(5) 


3.71(5) 


0.637(4) 


0.480(3) 
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5.271 


0.1412 


0.350(1) 


2.94(4) 


2.71(3) 


1.511(7) 


0.9282(3) 
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0.426(2) 


2.92(7) 


2.21(9) 


1.076(3) 


0.693(10) 
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0.344(2) 


3.27(28) 


2.50(8) 


0.933(15) 


0.566(1) 


10 
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0.0304 


0.347(2) 
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3.85(7) 


2.63(1) 


1.081(1) 


0.803(1) 
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5.58 


0.096 


0.2716(10) 


4.8(1) 


2.94(9) 


0.9513(8) 


0.777(2) 



Table 1: Different quantities along LCP2 (upper part) and the "non-LCP approach" (lower part). 
The middle part shows new simulation points that were used to determine LCP2. As usual, the 
numbers in the parenthesis indicate the estimated statistical errors of the last digits. The effective Nt 
values for the "non-LCP approach" were calculated as follows. We determined four different critical 
couplings on A''f=4,6,8 and 10 lattices, while the quark mass was fixed at ma=0.096. Effective Nt-s 
were determined by interpolating between and slightly extrapolating from these points. 

By the finite temperature technique, described above, only a few points of the LCP's 
can be obtained. To interpolate /3(a) between these points (and extrapolate slightly away 
from them) we use the renormalization group inspired ansatz proposed by Allton [47]. Note 
that not only the quark mass parameter (m^a) can be used as a parameter of the LCP's, 
but a particularly illustrative parametrization is obtained by inverting eq. ( |3.lD and using 
Nt as a continuous parameter. (We use later this Nt parametrization for constructing the 
LCP on the fi — P plane by a linear combination of LCPi and LCP2, see Section 5.) Figure 
I shows LCPi and LCP2 with our simulation points. A line of constant physics obtained 
by renormalized quantities (LCPg) and the simulation points in the "non-LCP" approach 
are also shown. Note that even though the determination of the LCPi and LCP2 are done 
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3.2 


5.19 


0.096 


1.197(3) 


2.53(1) 


0.908(3) 


0.529(3) 


4 


5.271 


0.096 


1.190(8) 


2.87(2) 


0.874(3) 


0.552(4) 


4.8 


5.36 


0.096 


1.18(6) 


3.35(18) 


0.835(6) 


0.612(3) 


6.21 


5.5 


0.096 


1.25(3) 


4.16(8) 


0.852(11) 


0.743(2) 
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5.58 


0.096 


1.30(3) 


4.57(10) 


0.80(3) 


0.817(3) 



Table 2: Diincnsionless combinations along LCP2 (upper part) and the "non-LCP approach" 
(lower part). The middle part shows new simulation points that were used to determine LCP2. 
Quark mass dependent dimensionless quantities change more in the "non-LCP approach" . 



on finite temperature lattices, the obtained bare parameters are used in the rest of the 
paper for T = and T ^ simulations. 

4. Equation of state along the lines of constant physics (LCP) at fj. ~ 

In previous studies of the EoS with staggered quark actions, the pressure and the energy 
density were determined as functions of the temperature for fixed value of the bare quark 
mass mqa in the lattice action [31-33]. In these studies a fixed Nt was used (e.g. Nt = 4: 
or 6) at different temperatures. Since T = l/{Nta) the temperature is set by the lattice 
spacing which changes with (3. This convenient, fixed bare m^a choice leads to a system 
which has larger and larger physical quark masses at decreasing lattice spacings (thus, at 
increasing temperatures). Increasing physical quark masses with increasing temperatures 
could result in systematic errors of the EoS. 

Clearly, instead of this sort of analysis (in the rest of the paper we refer to it as "non- 
LCP approach") one intends to study the temperature dependence of a system with fixed 
physical observables, therefore on an LCP. 

Recently, the EoS was also studied by using an action with dynamical Wilson quarks 
[48] . In this case there is no convenient choice similar to the approach with fixed niqa { "non- 
LCP approach"). The authors determined the lines of constant physics and the /? functions. 
They obtained the EoS along different LCP's in the range of mj^/mp = 0.65 — 0.95. 

In our analysis we use full QCD with staggered quarks along the LCP and compare 
these results with those of the "non-LCP approach" . 
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Figure 2: The lines of constant physics (LCPi and LCP2) on the (3 vs. mudd plane. The 
strange quark mass is given by = 2.08m„(i for both LCP's. The simulation points are shown 
by squares/triangles and connected by dashed/dotted lines for LCP1/LCP2, respectively. An ap- 
proximate LCP obtained by renormalized quantities (LCP2, the definition is m-^lmp k, 0.626) is 
shown by a solid line. The diamonds along a horizontal line represent the simulation points in the 
"non-LCP" approach. Additional 3 diamonds in the vertical direction show the simulation points 
used to test the path independence of the integral method (see next Section). 

In order to be self-contained in this Section we review the = technique to determine 
the interaction measure (e — 3p) by using the non-perturbative /^-function. Then we shortly 
discuss the integral method [34] to determine the pressure. In practice, the energy density, 
entropy density or the speed of sound can be directly deduced from the interaction measure 
and the pressure (c.f. sT = e + p and = dp/de). We review the basic formulas and 
emphasize the issues related to the EoS determination along an LCP. 

The energy density and pressure are defined in terms of the free-energy density (/): 



Expressing the free energy in terms of the partition function (/ = —T/V log Z = 
-Td{\ogZ)/dV) we have: 




p{T) 



(4.1) 



e(r) 



91og Z 
V dT 



p{T) = T 



dV 



(4.2) 
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As mentioned before, we use the same lattice spacings in the temporal and spatial 
directions. The temperature and volume are connected to this lattice spacing by 



aNt 



V = a-^Nt. 



(4.3) 



By varying the bare parameters one can change the lattice spacing; however, it is not 
possible to change T and V independently. This is the reason why eqs. (4!2|,4.3) cannot be 
used directly. Instead one determines the interaction measure by the help of the (3 function 
and the pressure by using the integral method. 

we see that the interaction measure (e — 3p)/T^ is directly 



Inspecting eqs. (|4.2| , 
proportional the total derivative of log Z with respect to the lattice spacing: 



3p 



Nf d{\ogZ) 

m'' da ■ 



(4.4) 



Here, the derivative with respect to a is defined along the LCP, which means that only 
the lattice spacing changes and the physics (in our case niq/Tc) remains the same. The 
variation of the lattice spacing can be controlled by the bare parameters of the action (/? 
and mqa), so we can write: 



d_ 

da 



'dadp 



d{mga) d 
da d{niga) 



(4.5) 



The partial derivatives with respect to a should be taken along the LCP. Since the LCP 
is defined by rUq/Tc = const., the partial derivative d{mqa) / da becomes simply ruq. The 
derivatives of log Z with respect to (3 and mq are the plaquette and ^^q averages multiplied 
by the lattice volume. We get: 



3p 



2^4 




LCP 



+ ^^qmq j 

q / 



(4.6) 



The pressure is the other basic quantity, which is usually determined by the integral 
method [34]. The pressure is simply proportional to logZ, however it cannot be measured 
directly. One can determine its partial derivatives with respect to the bare parameters. 
Thus, we can write: 



P_ 

2^4 



{f3,mqa) 



s i(/3o,r7iqoa) 



d{(3, mqO) 



dlogZ/dp 
9 log Z/d{mqa) 



rp4 • 



(4.7) 



Since the integrand is the gradient of log Z, the result is by definition independent of the 
integration path. We need the pressure along the LCP, thus it is convenient to measure 
the derivatives of log Z along the LCP and perform the integration over this line. For the 
subtracted vacuum term we used the zero temperature pressure, i.e. the same integral on 
Nfo = Ng lattices. The lower limits of the integrations (indicated by /?o and m^o) were set 
sufficiently below the transition point. By this choice the pressure becomes independent 
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of the starting point (in other words it vanishes at vanishing temperature). In the case of 
2 + 1 staggered QCD eq. ( [4.7D can be rewritten appropriately and the pressure is given by 

( ^^^^ \ 
i^^ud) , (4.8) 

where we use the following notation for subtracting the vacuum term: 



- -N^ / d{P,muda,rnsa) 



rp4 



{l3Q,mqoa) 



{0{f3, m)) = 0(13, m)y_,o - 0(/3, m)^^,. (4.9) 

The integral method was originally introduced for the pure gauge case for which the 
integral is one dimensional, it is performed along the (3 axis. Previous studies for staggered 
dynamical QCD [31-33] used a one-dimensional parameter space. Note that for full QCD 
the integration should be performed along a path in a multi-dimensional parameter space. 
In our 2 + 1 flavour staggered QCD case the parameter space is three-dimensional. 

According to eq.(4.7) the gradient of logZ is measured, thus the integral should be 



independent of the path. We explicitely checked this independence by performing the 
integration along different paths. For this purpose we determined the normalised pressure 
{jp/T^) at /? = 5.5, rriuii = 0.048 (this point corresponds to T = 2Tc with quark masses 
defined by eq.( |3.lD ). In the calculation two different integration paths were used. In the 
first case we integrated along the LCP. We obtained p/T^ = 6.15(5). In the second case 
the integration path contained two pieces (see the path defined by the diamonds of Figure 
|2|). The first piece is a /3 integral at constant m^, the second one is a rriuda 'tUs integral (with 
fixed mass ratio) at constant /3 value. For the second path we obtained p/T'^ = 6.15(4). 
As it can be seen the two different paths give the same result for the pressure (within the 
estimated uncertainties). 

The EoS is determined along the LCP's and also by using the "non-LCP approach". 
The results along LCP2 and LCPi are very close to each other (the T = pion masses differ 
only by 10% for the two LCP's; for the m-,^/mp dependence of the EoS with Wilson quarks 
see [48]). Figure ^ shows the EoS at vanishing chemical potential on A'^t = 4 lattices for 
LCP2 and for the non-LCP approach. The pressure and e — 2>p are presented as a function of 
the temperature. The parameters of LCP2 and those of the non-LCP approaches coincide 
at T = Tc- The Stefan-Boltzmann limit valid for Nt = 4: lattices is also shown. 

It is interesting to compare the pressure difference between LCP and non-LCP ap- 
proaches at finite temperature in case of the ideal and interacting gas. In the ideal gas 
case there is a relative 2.7% difference between the pressures, while this difference is three 
times larger in the interacting case (9.4%) at T = 2.5 • Tc- This means that considering 
only the mass-effect of the ideal gas, we would make a significant underestimation for the 
interacting case. 

5. Reliability of reweighting 

In order to determine the EoS for the 7^ case (when direct simulations are not possible) 
we use the multi-parameter reweighting proposed in [17] The aim of this section is to 
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Figure 3: The equation of state at fj,=0. (a) The left panel shows the pressure p, as a function 
of the temperature. All quantities are normalised by T*. The solid line connects the data points 
obtained along the LCP2, whereas the dashed line connects the data points obtained in the "non- 
LCP approach" (see text). The Stefan-Boltzmann limit is also shown by an arrow for Nt=4: lattices. 
The EoS along the LCP and the "non-LCP approach" differ from each other at high T. (b) The 
same for the "interaction measure" (e-3p). 

point out that the multi-parameter reweighting [17] is reliable. (For another study of 
reweighting see [49].) We also determine its region of validity by a suitably estimated error. 
Furthermore we show that in certain cases the multi-parameter reweighting is much more 
reliable than single-parameter techniques (Glasgow- type reweighting [18]) in the sense that 
one can reach farther regions. On top of this we demonstrate that the major requirement to 
be met by the sample used for the reweighting is that it should be as variable as possible. 
Variability here means that it should contain configurations from the various phases or 
configurations from the different Z(3) sectors in case of imaginary chemical potentials. To 
start with let us briefly review the multi-parameter reweighting. 

As proposed in [17] one can identically rewrite the partition function in the form: 

Z(m,^,/5) = jvUexp[-Sbosif3o,U)]detMimo,fi = 0,U) (5.1) 

/ \ <^ (R Tn^<^ (R TTM detM(m,/i,[/) ) 
^eM-SUP, U) + SUPo, ^)] detM(mo,/. = 0,[/) / ' 

where U denotes the gauge field links and M is the fermion matrix . The chemical 
potential n is included as exp(a/x) and exp(— a^) multiplicative factors of the forward and 
backward timelike links, respectively. In this approach we treat the terms in the curly 
bracket as an observable - which is measured on each independent configuration, and can 

^For n/ 7^ 4 staggered dynamical QCD one simply takes fractional powers of the fermion determinant. 
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be interpreted as a weight - and the rest as the measure. Thus the simulation can be 
performed at /i = and at some Pq and tuq values (Monte-Carlo parameter set). By 
using the reweighting formula (|5.1[) one obtains the partition function at another set of 
parameters, thus at /i / 0, /? 7^ /3o or even at m / niQ (target parameter set). 

Expectation values of observables can be determined by the above technique. In terms 
of the weights (i.e. the expression in the curly bracket of eq. (^.1|)) the averages can be 
determined as: 

^ E{^(/3^/|,m ^)}0(A/^,^,^) _ (5,2) 

It is clear that simulating at a given Monte Carlo parameter set the errors increase as 
we go farther and farther with the target parameter set. We need an error estimate of the 
reweighted quantities which shows how far we can reweight in the target parameter space. 
This way we could determine the borderline of the region outside of which the reweighting 
already fails to give reliable predictions. To achieve this we checked several error definitions 
including the jackknife and the statistical errors described in the papers [50, 51]. Both 
proved little to draw the limit of the reweighting procedure. The reason for this is that 
neither of them contains the systematic errors occuring because of the finite sample size. 
We note that in [51] there is an error estimate recommended for taking the finite sample 
size into consideration. However, this estimate uses an approximation of the distribution 
of the sample. An approximation of this kind can not be justified in case of staggered 
dynamical QCD. 

Eventually, with a new technique different from the ones mentioned above we succeeded 
to define a reliable error estimate for the reweighted quantities and to draw the limit of 
the reweighting procedure. After defining the new technique we show its application on an 
example. Then we demonstrate the problems that crop up while using the jackknife error 
and furthermore show how our new method overcomes these difficulties. 

The steps of the new procedure are as follows. First we assign to each and every 
starting configuration the weight w valid at the chosen target parameter set i.e. 



w 



exp |a/3 • y • {PI) + ^[lndetM(/i) - Indet Af(;U = 0)]} (5.3) 



where V is the lattice volume. Note that only the /3 and the /i parameters differ from 
their values at the simulation point. We generate a new set of configurations using these 
weights for Metropolis-like accept-reject steps. Let the original and new configuartions be 
<I>j and $^ and the corresponding weights (according to eq. (^.3|)) Wi and w'^. The first new 
configuration is = $1. Then for i > 1, <I>^ = <I>j with the probability max{l, } and 
$^ = otherwise. For complex weights we have to take either the absolute values or the 
absolute values of the real parts of the weights. The configurations of the new sample 
are taken with a unit weight to calculate expectation values, variances and integrated 
autocorrelation times [50] (the latter grows at every rejection due to the repetition of 
certain configurations). If our initial sample consisted of a large number of independent 
configurations then for real weights this method regenerates a sample which is theoretically 
perfect for the target parameters. For complex weights the new set of configurations can 
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not be used to measure observables, but it can be used for getting an error estimate. 
Therefore, in general, we will use eq. ( |5.2| ) to get the expectation values and the newly 
generated sample to determine the errors (see later). 

We still have to clarify two important points before our new procedure is finalised. 
The first is the question of thermalization. The second question is (in an extreme case): 
what happens if there is no configuration in the starting sample which would be relevant at 
the target parameters. We solved the first problem by defining a thermalization segment 
at the beginning of every newly generated sample which we cut off from the sample before 
calculating the expectation values and the errors. An obvious assumption is to claim 
that the already thermalized sample containing valuable information starts with the first 
different configuration right after the one with the largest weight. This ensures that if there 
is only one configuration in the initial sample which "counts" at the target parameters then 
this information will not be lost. The second problem can not be solved perfectly. The 
problem occurs e.g. when a phase transition is very strong, that is the physical quantities 
change significantly during the transition (first order transition) and we would like to 
reweight starting from deep inside in one phase into deep inside into the other phase. Then 
the solution can only be a huge sample size (or small lattice size) which allows the presence 
of a few configurations of the target phase in the initial sample. It is, however, difficult to 
define what a huge sample really means in a general, practical way. 

To illustrate the new technique let us take a look at the nj = 4 flavour case at 
niqa = 0.05 bare quark mass on 4-6^ size lattice at imaginary chemical potential. Note that 
for purely imaginary chemical potentials direct simulation is possible, therefore it is possible 
to check the validity of any error estimation method. Direct simulation is thus used as a 
standard. Accordingly, we determine the plaquette {PI) and the ipip expectation values and 
their uncertainties at /5 = 5.085 and Im(/i) ^ 0, that is at the target, imaginary chemical 
potential values using three different techniques. First technique is direct simulation (60000 
configurations). The other two methods are based on reweighting Im(/i) = data. We 
carried out simulations at Im(/x) = in the phase transition point, i.e. at /3 = 5.04 (~ 2500 
independent configurations) and at /3 = 5.085 (~ 7000 independent configurations)^. From 
these two starting points with the use of the reweighting we predicted the plaquette {PI) 
and the ^jJ^ expectation values and their uncertainties at /3 = 5.085 and Im(/x) ^ 0, that is 
at the target. In the following we refer to the technique starting from the phase transition 
point at Im(/i) = as multi-parameter reweighting, while the technique starting from 
P = 5.085 (fixed at the target value) and Im(/x) = as Glasgow-type reweighting. We used 
the new Metropolis-type method defined above in both cases. As a check for both starting 
points we also calculated the plaquette and the -^V expectation values directly by ( |5.2| ) at 
the target points, and observed that the results are very similar to the results obtained by 
the Metropolis-type method. 

The expectation values and their errors obtained by the new (Metropolis-type) method 
along with the direct simulation points are shown on the left panel of Figure |^. When 
making the figure we had to make use of an additional information, namely the fact that 

*These parameters are identical to the ones used in [17] 
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Figure 4: (a) The left panel shows '0'!/' expectation values and error estimates of three dif- 
ferent methods. The squares correspond to direct simulations (out of 60000 configurations). The 
crosses denote the results of the multi-parameter reweighting method (out of « 1200 independent 
configurations) while the circles are the points of the Glasgow-type reweighting (out of « 1200 
independent configurations). In case of the latter two methods the error estimates are determined 
using the new Metropolis-type method, (b) In the right panel the meaning of the symbols are 
unchanged but the sample sizes for the reweighting techniques are increased to 2500 independent 
configurations in case of the multi-parameter reweighting and to 7000 in case of the Glasgow-type 
reweighting. In case of the multi-parameter reweighting the increased sample size only decreases 
the error estimate, maintaining consistency with the direct simulation result. On the other hand 
in case of the Glasgow-type reweighting the fact that the small sample (used in the left panel) does 
not contain any configuration from the target phase causes systematic errors in the expectation 
values and mainly in their uncertainties. This is because in case of strong phase transitions ~ like 
in our case - it is far too difficult to define a reliable thermalization stage by virtue of the weights. 
Increasing the sample size helps in this situation leading to the results of the right panel. The very 
large error estimates in the Glasgow case are realistic, the prediction does agree within errors with 
the direct simulation result. 

the /3— Im(/i) plane consists of three sectors. If we want to make predictions about all of 
them based on the data at = we can do it by explicitely using the Z(3) symmetry. That 
is we use the starting sample in such a way that every plaquette value is included three 
times: once with the original weight and also with the weights of shifted to fi^in/Nt/S 
arguments.^ Performing the reweighting like this we obtained Figure ^. In the left panel 
the second problem mentioned in the previous paragraph is seen. I.e. in case of Glasgow- 
type reweighting the small sample (1200 independent configurations) causes the problem. 
Namely, it does not contain any configuration belonging to the target phase when the target 
chemical potential is around Im(//) = 7r/12. Then even our new method gives misleading 
results. Increasing the samples size, as soon as a single, dominant configuration turns up, 

^Note that Z(3) symmetry is a specific feature of QCD at imaginary chemical potential, it has nothing 
to do with our new Metropolis-type method. 
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the new method supphes us with rehable error estimates. (See the right panel of Figure 
where we used 7000 independent configurations in place of the 1200 configurations that 
lead to the wrong results on the left panel.) The right panel truly reflects the correctness of 
the multi-parameter reweighting in contrast to the single-parameter reweighting Glasgow- 
method. The infinitely large errors of the Glasgow-method indicate that the whole sample 
is thermalization, that is it does not provide information about the expectation value in 
the corresponding imaginary chemical potential points. 
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Figure 5: Plaquette expectation values obtained by direct simulations and reweighting at various 
/3-s. The squares and their errorbars denote the results of direct simulation. The initial sample used 
for reweighting consisted of 33000 configurations. The parameters were: f3 = 5.274, the critical f3 
in the 71/ = 2 + 1 flavour case when mu.d = 0.096, = 2.08mu.d are the quark masses (in 
lattice units) on a 4 • 8"^ size lattice. The circles and the errorbars belonging to them represent 
the newly introduced Metropolis-type method. Infinitely large errors mean that the whole, newly 
generated sample is thermalization so it does not contain information about the distribution at 
the (3 parameter in question. Tlic crosses refer to the values given by the formula (5.2), the errors 
assigned to them are jackknife errors. It is clear that the jackknife errors do not notice the limit of 
validity of the reweighting procedure, but for small |A/3|-s they provide a good error estimate. 

An example of the troubles occuring when using the jackknife error estimate will be 
shown in case of single-parameter reweighting in the (3 parameter. Thus one can use the 
new (expensive) error estimate to provide the limit of the applicability of the reweighting 
procedure after which the simpler jackknife method can be applied in the appropriate 
region. (This strategy was followed in the rest of the paper to calculate the EoS.) In our 
example we take 4 • 8^ size lattice in case of ?7i„ ^ = 0.096, = 2.08mu^d bare quark 
masses (in lattice units) at the critical point, i.e. at (3 = 5.274. A sample simulated in 
this point consisting of 33000 configurations was reweighted by (|5.2| ). We evaluated the 
errors of the reweighted quantities by the Metropolis-type and by the jackknife methods. 
These are shown together with the reweighted plaquette expectation values in Figure |5|. 
It can be seen that the new method gives information on the errors and also on the limit 
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of the applicability of the reweighting procedure. On the other hand the jackknife error 
can not be used outside a restricted region and does not provide any clue on the region of 
applicability of the reweighting procedure. 

What can we say about the errors in the real fi case based on the new error estimation 
procedure? The Metropolis steps can only be taken using either the absolute values or 
the absolute values of the real parts of the weights. It is clear that both possibilities 
represent only approximations to the true reweighting given by eq. ( |5.2| ). Whereas the 
expectation values provided by equation (|5.2| ) are exact in principle, we may test the above 
defined Metropolis-type procedures for the errors provided by them. We can claim based on 
b igure I that on a 4 • 8^ size lattice in the fi £ [0, 0.3] region and also in all other examined 
cases the three methods do not provide significantly different expectation values, differences 
are tolerated in the errors. So we assign the errors of the Metropolis-type method given 
through the use of the absolute values of the real parts of the weights to the expectation 
values obtained from equation ( [5.21) . 
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Figure 6: Average plaquette versus the real chemical potential. 33000 configurations were 
simulated at the parameter set: (3 = 5.274, m„^d = 0.096, nis = 2.08mu.d on a 4 • 8^ size lattice. 
This is at the critical f3 in the Uf = 2 + 1 flavour case. The circles refer to the plaquette values 
given by (5.2), their errorbars indicate jackknife errors. The squares belong to a Metropolis-type 
method using the absolute values of the real part of the weights while in case of the triangles the 
absolute values of the weights were used. In both cases the errors are given by the square roots of 
the variances times the autocorrelation time. 



6. Best reweighting lines and the LCP's 

We can define reweighting lines on the (3-^ plane so as to make the least possible mistake 
during the reweighting procedure. To do this we introduce the notion of overlap measure 
which we denote by a. The overlap measure is the normalised number of different config- 
urations in the sample created with the Metropolis-type reweighting after cutting off the 
thermalization. We plotted the contour lines of a in the left panel of Figure 0. The dotted 
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areas are unattainable, that means here the overlaps vanish, the errors are infinitely large. 
The best reweighting lines can be defined for each simulation point. For a given value of 
we choose /? corresponding to the maximal value of a. The points of the best reweighting 
line are given by the rightmost points of the contours of constant overlap in Figure ^ a. 




Figure 7: (a) The left panel shows the real chemical potential-/? plane. 33000 configurations 
were simulated at the parameter set: /3 = 5.274, rriu^d = 0.096, mj — 2.08m„_d on a 4 • 8'^ size 
lattice. This is at the critical /3 in the n/ = 2+1 flavour case. The dotted lines are the contours 
of the constant overlap. The dotted area is the unknown territory where the overlap vanishes, i.e. 
the sample regenerated from the initial one consists of nothing but thermalization. The roughness 
of the border separating the unknown territory is due to the presence of random numbers in the 
procedure. The solid line is the phase transition line determined by the peaks of susceptibility, (b) 
In the right panel the volume and the /i dependence of the overlap (a) is shown. Upper curves 
correspond to smaller lattice sizes, 4-6^, 4-8'^, 4 • 10"^ and 4 • 12^ respectively. The half width (/Lti/2; 
defined by a |^^^^ = 1/2) scales according to: ^1/2 oc with 7 w 1/3. 

It is important to examine the volume dependence of reweighting by using the notion 
of the overlap. The right panel of Figure |^ shows the /i dependence of the overlap at fixed 
(3 and quark mass parameters for different volumes (V^ = 4 • 6^, 4-8^, 4 • 10'^ and 4 • 12'^). 
As expected, for fixed larger volumes result in worse overlap. One can define the "half- 
width" {^112) of the dependence by the chemical potential value at which a = 1/2. One 
observes an approximate scaling behaviour for the half-width: /Ui/2 oc V'"^ with 7 ~ 1/3, 
which is much better than a crude XjyfV estimate. Thus the volume dependence of — /? 
reweighting is less restrictive than that of reweightings in other parameters. 

Expectation values of observables at finite chemical potential are determined according 
to eq. ( p.2p . There are quantities (e.g. plaquette, Wilson-loop, Polyakov-line) for which the 
\x reweighting means only a change in the weight of the configuration but not in the value 
of the observable. For other quantites (e.g. chiral condensate, fermion number) not only 
the weight but also the value is infiuenced by the chemical potential since they are directly 
expressed by the fermion matrix. It is instructive to study the volume dependence of both 
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Figure 8: (a) The /i dependence of the plaquette (left panel). Results obtained on 4 • 8'^ lattices 
are shown by horizontally, those on 4 • 10^ lattices are shown by vertically and those on 4 • 12^ 
are shown by diagonally lined bands. The results agree within their statistical errors, (b) The fi 
dependence of fermion number density (right panel). The difference of 4-8^, 4-10^ and 4-12'^ results 
cannot be resolved (solid line). The dashed line shows the Stefan-Boltzman limit. Simulations were 
done at /3 = 5.35. 

types of these quantities. The left panel of Figure |8| shows the average plaquette, the right 
panel the fermion number density, as functions of the chemical potential. Both quantites 
are determined on 4 • 8^, 4 • 10^ and 4 • 12^ lattices. The volume dependence is practically 
negligible for the fermion number density, whereas it is moderate for the average plaquette. 

It is obvious that the two-parameter reweighting used in the previous section does 
not follow the LCP (/3 gets smaller but the quark mass remains moa). We show in the 
following that the best reweighting line along the LCP can be practically determined by 
two techniques: a. three-parameter reweighting, b. interpolating method. 

a. As it can be seen in the left panel of Figure |^ the change in /? is not very large for 
the two-parameter reweighting. Therefore, one can remain on the LCP by a simultaneous, 
small change of the mass parameter of the lattice action. This results in a three-parameter 
reweighting (reweighting in ma, j3 and no). For some fixed fJ-target two constraints are 
needed to determine mtargeto- and Ptarget- One of them is the generalization of the new. 
Metropolis-type reweighting to some nritargeto- 7^ rriQa and the other one is the LCP, which 
relates^'' the mass parameter to the gauge coupling, i.e. one has to fulfill the following 
condition 

(^target = PLCP{mtargeta) ■ (6.1) 

Similarly to the two-parameter reweighting one can construct the best three-parameter 
reweighting line, which is shown in Figure ^. 

^"in our specific case the PLCpima) relationship is given in Figure El 
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Figure 9: Different types of best weight lines are shown on the /i-/3 plane. Three-parameter 
(m, (3 and fi) reweighting is used for the parameter set of LCP2(A^t — 4) (dashed line). Along the 
line the system remains on the LCP, the small change in (3 is compensated by the change of the 
mass parameter. Two parameter reweighting is used for the parameter set of LCPi(iV( = 4) (upper 
solid line) and for that of LCP2{Nt — 4) (lower solid line). The interpolating technique is used 
to determine the line of constant physics on the plane (dotted line). Note that the results of 
the three-parameter reweighting technique and that of the interpolating technique agree within the 
statistical uncertainties. 

b. The other possibility to stay on the line of constant physics at finite fi is the interpo- 
lating technique. One uses the two-parameter reweighting for two LCP's and interpolates 
between them. As we discussed in Section 3 not only the mass parameter, but also Nt (as 
a continuous parameter) can be used to parametrize the LCP's. Thus the data of Figure 
I give (5i{Nt,n = 0) for LCPi and l32{Nt,fi = 0) for LCP2, respectively. The same Nt 
values for the two LCP's correspond to different mass parameters, which can be written 
as (am)i = {am)i{f3i, Nt, fj, = 0) and (am)2 = {am)2{(32, Nt, fj, = 0). Note that the two- 
parameter reweighting to finite /x does not change the mass parameter. Let us take some 
fixed Nt value and perform two-parameter reweightings for both /3i and f32 parameter sets. 
This results in /3i(/x) and /32(a') best weight lines (these functions for Nt=^ are shown in 
Figure]^ by solid lines). Though these individual curves leave the LCP-s, there is a way to 
stay on e.g. LCP2 by interpolating between /3i(^) and /32(/^)- Since the mass parameters of 
LCPi and LCP2 are close to each other we use a linear interpolation between them. The 
obtained gauge coupling is /3i(/u) and the mass parameter is {am)i{ij,). 

A(/i) = C/3i(^) + (l-C)/32(/i), (6.2) 
(am)i(/i) = Ciam)i{(3i,Nt,fi = 0) + (1 - C)(am)2(/52, iV*, = 0), 

A = (32{{am)i). 

The first two equations define the interpolation and the third condition guarantees that we 
stay on LCP2. These three equations determine the three unknown variables: Pi, {am)i 
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and the interpolating parameter The best weight hne obtained by the interpolation 
technique is also shown in Figure As it can be seen the result of this method and the 
predicition of the three-parameter reweighting agree quite well. This indicates that the 
requirement for the best overlap selects the same weight lines even for rather different 
methods. 



7. Equation of state at non-vanishing chemical potential 

In this section we study the EoS at finite chemical potential. Since we are interested in 
the physics of finite baryon density we use = = /^d / for the two light quarks and 
fis = for the strange quark. The same technique and the same equations could be used 
for the fj,s ^ case (one should simply add the strange quark contribution to that of the 
light quarks). 

The energy density and pressure can be derived in a similar way as in the case of 
vanishing chemical potential discussed in Section ^. However, now the grand canonical 
potential (lo) is used instead of the free energy. The definitions are: 



UJ 



T— 

dT 



duj ,^ . 



(7.1) 



Expressing the grand canonical potential in terms of the grand canonical partition 
function (oj = —T/V log Zgc) we get similar results as in the = case. In order to 
make this transparent we explain the calculation in detail. As usual one has to separate 
the space-like (a^) and time-like {at) lattice parameters. This way the derivatives with 
respect to T, V and fi can be replaced by derivatives with respect to as, and fx, where 
^ = at/as, fi = nag. It does not matter whether we write fiat or /uog, because we will set 
^ = 1 at the end of calculation. The new (lattice) variables can be expressed by the old 
(thermodynamical) ones: 
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And the derivatives: 



Ns 1 



1 ..yi/3. 



_d_ 

dT 

d_ 
dV 

d_ 



TV 



Us, A* 




0.8, i 



In particular, for e — 3p we have: 
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Now we replace thermodynamical variables to lattice variables 
e — 3p 



T d log Zg, 



T 
V 



dlogZ, 



fiT dlog Zgc 
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dlogZ, 
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Setting ^ = 1 we get: 



e-3p _ 9(log Z, 



gc) 



(7.3) 



(7.4) 



^a=const 

We can again write the derivative with respect to a in terms of derivatives with respect 
to the bare parameters. The additional bare parameter a/x is kept constant, so it will not 
generate an extra term: 
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(7.5) 



IJ,a=const 



We have the same expression for e — 3p as in the fi = case, the only difference is that 
now all observables should be evaluated at finite fi: 
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For the pressure we have an additional variable in the integral, the gradient of log Zgc 
has now an extra component in the afi direction: 
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The subtracted po term is now the pressure at T = and fJ, = 0, i.e. the same as before. 
We can rewrite the above equation in terms of observables: 
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where the superscript /i denotes the expectation values of the observables at reweighted p, 
values: 

{0{P, /i, m))('^) = 0{P, /i, m)^^o - 0{P, p = 0, m)^^^. (7.9) 

The result should again be independent of the integration path. We chose the following 
paths (cf. Fig. IC). First we integrated from (/3o, aniio) to some (/3i, anin) along the LCP 



- 21 - 




Figure 10: Illustration of the integral method at finite chemical potential (left panel). The solid 
lines are /i = const lines on the fia — P plane. Dashed lines are the best reweighting lines starting 
from different simulation points. Arrows show the path of integration we used when evaluating 
eq. (7.7). 



with fj, = 0. Then we followed the line of best reweighting to reach the required {(3, ami, 
ajj,) point. We checked that the result remains unchanged if we do it in the opposite way, 
i.e. we go first along the best reweighting line and then along the constant a// line. The 
results along the two integration lines were equal within statistical errors. 

We present lattice results on Ap{fi,T) = p{fi / 0, T) — p{fj. = 0, T), e{fj,,T)-3p{fj,,T) 
and UBifJ-jT). Our statistical errorbars are also shown. They are rather small, in many 
cases they are even smaller than the thickness of the lines. 

On the left panel of Fig. |l^ we present Ap/T'^ for five different fj, values. On the 
right panel normalisation is done by Ap'^^, which is Ap{fi, T — > oo). Notice the interesting 
scaling behaviour. Ap/Ap'^^ depends only on T and it is practically independent from 
in the analysed region. The left panel of Fig. 12 shows the same pressure difference as 
a function of ^b/T for five different temperatures. The right panel shows dimensionless 
baryon number density (ub/T^) as a function of fJ-s/T at the same temperatures. The left 
panel of Fig. |l^ shows e-Sp normalised by T^, which tends to zero for large T. The right 
panel of Fig. 13 gives the dimensionless baryonic density as a function of T/Tc for different 
H-s. 

Here we summarise what kind of errors occur during our lattice calculations. The error 
coming from reweighting has been amply discussed in Section 5. Another source of error 
is the finiteness of the physical volume. The volume dependence of physical observables 
is smaller than the statistical errors for the plaquette average or quark number density. 
Systematic errors coming from non-zero microcanonical step-size cause a 0.2% and 0.1% 
error in the value of physical observables for zero temperature and finite temperature 
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Figure 11: The equation of state at, = 100, 210, 330, 410 and 530 MeV. The left panel shows 
the pressure difference between the fi — and /i ^ cases normalized by T**, whereas on the right 
panel the normalization is done by the A^f = 4 lattice Stefan-Boltzman limit (in continuum this 
limit is n/(/i/T)^). Note that Ap/Ap^^ seems to show some scaling behaviour (it depends mostly 
on the temperature; whereas its dependence on /i is much weaker). Thus, the fi dependence of Ap 
is almost completely given by the /i dependence of the free gas. 



calculations, respectively. 
8. Conclusions, outlook 

In this paper we studied the thermodynamical properties of QCD at finite chemical po- 
tential /i. We used the overlap enhancing multi-parameter reweighting method proposed 
by two of us [17] and its generalization to 2 -|- 1 flavour staggered QCD [21]. Our primary 
goal was to determine the equation of state (EoS) on the line of constant physics (LCP) 
at finite temperature and chemical potential. 

We have pointed out that even at ij,=0 the EoS depends on the fact whether we are 
on an LCP or not. Note that previous results in the staggered formalism usually used the 
"non-LCP approach" . According to our findings pressure and e - 3p (interacion measure) 
on the LCP has different high temperature behaviour than in the "non-LCP approach". 

We discussed the reliabihty of the reweighting technique. We introduced an error 
estimate, which successfully shows the limits of the method yielding infinite errors in the 
parameter regions, where reweighting gives wrong results. We showed how to define and 
determine the best weight lines on the plane. 

We discussed the two-parameter and three-parameter reweighting techniques. Two 
techniques were presented (three-parameter reweighting and the interpolating method) to 
stay on the LCP even when reweighting to non-vanishing chemical potentials. 
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Figure 12: (a) The difference of the pressure at finite /is and /is = (left panel) as a function of 
Hb/T at T/Tc — 0.9, 0.98, 1.03, 1.2, 2.28. (b) Dimensionless baryon number density as a function 
of at T/Tc = 0.9, 0.98, 1.03, 1.2, 2.48. Dashed line shows the SB limit in both cases. The higher 
the temperature is, the closer the lines run to the SB limit. 
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Figure 13: (a) (e — 3p)/T^ (b) Dimensionless baryon number density as a function of T/Tc at 
Aii3=100,210,330,410 MeV and ^is^SSO MeV. 



We calculated the thermodynamic equations for ^ 7^ and determined the EoS along 
an LCP. We presented lattice data on the pressure, the interaction-measure and the baryon 
number density as a function of temperature and chemical potential. The physical range of 
our analysis extended upto 500 — 600 MeV in temperature and baryon chemical potential 
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as well. 

Clearly much more work is needed to get the final form of non-perturbative EoS of 
QCD. One has to extrapolate to zero step-size in the R-algorithm. Extrapolation to the 
thermodynamic and continuum limits is a very CPU demanding task in the 7^ case. 
Physical rriT^/mp ratio should be reached by decreasing the light quark mass. Finally, 
renormalised LCP's (LCP* in our notation) should be used when evaluating thermody- 
namic quantites. 
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